##Anderson Harbridge MIP Replication File.R
##Analysis of ideas from disproportionate info processing and motivated reasoning

##Budget data with pct change in each subaccount-year
load("Anderson_Harbridge_APR_Budget_Data.RData")
ls()

names(pct.change2)

#####################################################################################
##Data set of yearly percent changes in budget authority at the subaccount level.  ##
##Includes discretionary spending, including defense.                              ##
##Key variables are defined as follows:                                            ##
##nsa: subaccount identifier                                                       ##
##pctchange: real percentage change in budget authority from previous year         ##
##pctchangenom: nominal percentage change in budget authority from previous year   ##
##year: calendar year (1955-2002)                                                  ##
##genpurcat: Budget general purpose category (see Cogan 2002)                      ##
##secputcat: Budget secondary purpose category (see Cogan 2002)                    ##
##dayspastfy: Number of days past start of fiscal year that appropriations bill    ##
##   passed (NSA are nested within appropriations bills                            ##
##unemploy: Yearly unemployment rate                                               ##
##total_surplus_pct_lag: size of the surplus as a percentage of gross domestic     ##
##   product (GDP; lagged by 1 year)                                               ##
##paygo: 1 if a year when Pay-as-you-go was in place                               ##
##grh: 1 if a year when Graham-Rudman-Hollings was in place                        ##
##second.session: 1 if second session of Congress                                  ##
##net.dem: average of net seats gained by Democrats since previous year in House   ##
##   and Senate (as a percent)                                                     ##
##rep.senate: 1 if Senate controlled by Republicans                                ##
##rep.house: 1 if House controlled by Republicans                                  ##
##rep.pres: 1 if President is a Republican                                         ##
##s.no.dem: Number of Senate seats held by Democrats                               ##
##h.no.dem: Number of House seats held by Democrats                                ##
##s.no.rep: Number of Senate seats held by Republicans                             ##
##h.no.rep: Number of House seats held by Republicans                              ##
##new.event: 1 if new subaccount (inception)                                       ##
##zeroed.event: 1 if subaccount zeroed out                                         ##
##neg.change.event: 1 if negative percentage change of more than 3% (nominal)      ##
##pos.change.event: 1 if positive percentage change of more than 3% (nominal)      ##
##mess.with: 1 if percentage change is less than or equal +/-3% in magnitude       ##
##   (nominal) or if new subaccount, 0 otherwise                                   ##
##cut: 1 if (among subaccounts that are changed), neg.change.event=1 or            ##
##  zeroed.event=1, 0 if pos.change.event=1 or new.event=1, and NA if not changed  ##
##big.cut2: 1 if (among subaccounts that are cut), new.event=1 or percentage change##
##  is greater than 50% in magnitude (nominal), 0 if small cut, NA if increase or  ##
##  not changed (e.g. mess.with=0)                                                 ##
##big.increase2: 1 if (among subaccounts that are increased), subaccount is zeroed ##
##  out (zeroed.event=1) or percentage change is greater than 50% (nominal), 0 if  ##
##  small increase, NA if cut or not change (e.g., mess.with=0)                    ##
#####################################################################################


#############################################
##  Create new variables                   ##
#############################################
##regressions
##continuous variable
pct.change2$indicator2001<- ifelse(pct.change2$year==2001, 1, 0)

##Number of Democratic branches
pct.change2$no.dem.branch<- ifelse((pct.change2$rep.house==1 & pct.change2$rep.senate==1 & pct.change2$rep.pres==0) |
                                      (pct.change2$rep.house==0 & pct.change2$rep.senate==1 & pct.change2$rep.pres==1) |
                                      (pct.change2$rep.house==1 & pct.change2$rep.senate==0 & pct.change2$rep.pres==1), 1,
                                      ifelse((pct.change2$rep.house==1 & pct.change2$rep.senate==0 & pct.change2$rep.pres==0) |
                                      (pct.change2$rep.house==0 & pct.change2$rep.senate==1 & pct.change2$rep.pres==0) |
                                      (pct.change2$rep.house==0 & pct.change2$rep.senate==0 & pct.change2$rep.pres==1), 2,
                                      ifelse((pct.change2$rep.house==0 & pct.change2$rep.senate==0 & pct.change2$rep.pres==0), 3, NA)))

##Democratic president
pct.change2$dem.pres<- ifelse(pct.change2$rep.pres==0, 1, 0)

##Democratic Owned Issues
#4.4 housing assistance
#5.1 education
#5.2 job training
#6.4 housing
#6.5 environment
#6.8 energy and power subsidies
#6.9 community development aid
#8.7 mgmt of public domain

##Dem most owned
pct.change2$most.owned<- ifelse((pct.change2$genpurcat==4 & pct.change2$secpurcat==4) |
                                (pct.change2$genpurcat==5 & pct.change2$secpurcat==1) |
                                (pct.change2$genpurcat==5 & pct.change2$secpurcat==2) |
                                (pct.change2$genpurcat==6 & pct.change2$secpurcat==4) |
                                (pct.change2$genpurcat==6 & pct.change2$secpurcat==5) |
                                (pct.change2$genpurcat==6 & pct.change2$secpurcat==8) |
                                (pct.change2$genpurcat==6 & pct.change2$secpurcat==9) |
                                (pct.change2$genpurcat==8 & pct.change2$secpurcat==7), 1, 0)

##Add #dem chambers (H/S)
pct.change2$no.dem.chamber<- ifelse((pct.change2$rep.house==1 & pct.change2$rep.senate==1), 0,
                                    ifelse(((pct.change2$rep.house==1 & pct.change2$rep.senate==0) |
                                           (pct.change2$rep.house==0 & pct.change2$rep.senate==1)), 1,
                                           ifelse(pct.change2$rep.house==0 & pct.change2$rep.senate==0, 2, NA)))

##Average percent of seats held in House and Senate by Democrats
pct.change2$pct.dem<- (((pct.change2$h.no.dem/(pct.change2$h.no.dem + pct.change2$h.no.rep)) +
                           (pct.change2$s.no.dem/(pct.change2$s.no.dem + pct.change2$s.no.rep)))/2)*100
summary(pct.change2$pct.dem)

##separate for H and S
pct.change2$pct.dem.h<- (pct.change2$h.no.dem/(pct.change2$h.no.dem + pct.change2$h.no.rep))*100
summary(pct.change2$pct.dem.h)

pct.change2$pct.dem.s<- (pct.change2$s.no.dem/(pct.change2$s.no.dem + pct.change2$s.no.rep))*100
summary(pct.change2$pct.dem.s)



####################
##  Yearly Data   ##
####################
##load yearly data
yearly<- read.csv("yearly data.csv", header=T)

##define variables of interest
yearly$paygo<- ifelse(yearly$year >= 1991, 1, 0)
yearly$indicator2001<- ifelse(yearly$year==2001, 1, 0)

##Net seats gained by Democrats (as percent)
source('ar1.R')
yearly$net.dem <- ((((yearly$h.no.dem - mylag(yearly$h.no.dem, lag=1))/435)*100) + (((yearly$s.no.dem - mylag(yearly$s.no.dem, lag=1))/100)*100))/2
summary(yearly$net.dem)

##second session varaible
yearly$second.session<- ifelse(yearly$year %in% seq(1956,2002, by=2), 1, 0)

##number of Democratic branches
yearly$no.dem.branch<- ifelse((yearly$rep.house==1 & yearly$rep.senate==1 & yearly$rep.pres==0) |
                                       (yearly$rep.house==0 & yearly$rep.senate==1 & yearly$rep.pres==1) |
                                       (yearly$rep.house==1 & yearly$rep.senate==0 & yearly$rep.pres==1), 1,
                                       ifelse((yearly$rep.house==1 & yearly$rep.senate==0 & yearly$rep.pres==0) |
                                       (yearly$rep.house==0 & yearly$rep.senate==1 & yearly$rep.pres==0) |
                                       (yearly$rep.house==0 & yearly$rep.senate==0 & yearly$rep.pres==1), 2,
                                       ifelse((yearly$rep.house==0 & yearly$rep.senate==0 & yearly$rep.pres==0), 3, NA)))

##Democratic president
yearly$dem.pres<- ifelse(yearly$rep.pres==0, 1, 0)

##Add #dem chambers (H/S)
yearly$no.dem.chamber<- ifelse((yearly$rep.house==1 & yearly$rep.senate==1), 0,
                                    ifelse(((yearly$rep.house==1 & yearly$rep.senate==0) |
                                           (yearly$rep.house==0 & yearly$rep.senate==1)), 1,
                                           ifelse(yearly$rep.house==0 & yearly$rep.senate==0, 2, NA)))
summary(yearly$no.dem.chamber)

##Average percent of seats held by Democrats
yearly$pct.dem<- (((yearly$h.no.dem/(yearly$h.no.dem + yearly$h.no.rep)) +
                           (yearly$s.no.dem/(yearly$s.no.dem + yearly$s.no.rep)))/2)*100
summary(yearly$pct.dem)

##Gramm-Rudman-Hollings 1985-1990
yearly$grh<- ifelse(yearly$year >=1985 & yearly$year <= 1990, 1, 0)

#####################################################
##Most owned aggregated by year
yearly$mess.with.pct.mo<- NULL
i.index<- c(1955:2002)
for (i in 1:length(i.index)){
  yearly$mess.with.pct.mo[i]<- (sum(pct.change2$mess.with[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$mess.with[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$mess.with.pct.mo)

yearly$cut.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$cut.pct.mo[i]<- (sum(pct.change2$cut[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$cut[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$cut.pct.mo)

##for 50% as big
yearly$big.cut2.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.cut2.pct.mo[i]<- (sum(pct.change2$big.cut2[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.cut2[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.cut2.pct.mo)

##50% as big
yearly$big.increase2.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.increase2.pct.mo[i]<- (sum(pct.change2$big.increase2[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.increase2[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.increase2.pct.mo)


###########################################################
##           Distribution of cuts/increases              ##
###########################################################
##New analysis
##Distribution of small cuts, big cuts, small inc, big inc
##define small cuts and small inc
pct.change2$small.cut2<- ifelse(pct.change2$cut==1 & pct.change2$big.cut2==0, 1, 0)
pct.change2$small.increase2<- ifelse(pct.change2$cut==0 & pct.change2$big.increase2==0, 1, 0)

#####################################################
## Density of first and second changes for most owned
par(mfrow=c(1,1), las=1)
plot(density(pct.change2$pctchangenom[pct.change2$second.session==0 & pct.change2$most.owned==1], na.rm=T, from=-100, to=100), ##was -50 to 50
     ylim=c(0,.04),
     main="", #Distribution of Policy Changes
     xlab="Percent Change",
     ylab="Density")
lines(density(pct.change2$pctchangenom[pct.change2$second.session==1 & pct.change2$most.owned==1], na.rm=T, from=-100, to=100),
     lty=2,
     add=T)
legend(18, .04, c("First Session", "Second Session"), lty=c(1,2), bty="n")

############################
##Distribution of cuts
##Percent of cuts that are big
##Distribution of changes for Number Dem Branches

##Dem 1 branch
dem1.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==1 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==1 & pct.change2$most.owned==1], na.rm=T)))*100
dem1.pctbig.mo

dem2.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==2 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==2 & pct.change2$most.owned==1], na.rm=T)))*100
dem2.pctbig.mo

dem3.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==3 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==3 & pct.change2$most.owned==1], na.rm=T)))*100
dem3.pctbig.mo


##barplot
par(mfrow=c(1,1), las=1)
barplot(c(dem1.pctbig.mo, dem2.pctbig.mo, dem3.pctbig.mo),
         ylim=c(0,70),
     main="", #Distribution of Budgetary Cuts
     xlab="Number of Democratic Institutions",
        ylab="Percent of Big Cuts",
     names.arg=c("1", "2", "3"),
        cex.names=1)




######################################################
##        Analysis                                  ##
######################################################

############################################################
##Multi-level model with intercept varying by subaccount  ##
############################################################
library('Matrix')
library('lme4')


####################################
##With all subaccounts            ##
##interaction with second session ##
##Table 1                         ##
####################################
##Leave alone or change it
mod.no.dem2b.lmer.mess.with<- glmer(mess.with ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.mess.with)

##If change it, cut or increase it (1 is cut) (Note that table reverses to increase)
mod.no.dem2b.lmer.cut<- glmer(cut ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.cut)

##If cut, small % cut or big cut/zero out
##50% definition
mod.no.dem2b.lmer.big.cut2<- glmer(big.cut2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.cut2)

##If increase, small % or big/inception
##50% definition
mod.no.dem2b.lmer.big.increase2<- glmer(big.increase2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.increase2)

###############################################
##Use forge version to rpedict move from 1 Dem branch to 3 Dem branches
##copy complete data for regression
##Cut
foo.cut<- data.frame(cbind(pct.change2$cut, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.cut)<- c("cut", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.cut<- na.omit(foo.cut)

###############
##If all moved to first session
foo.complete.cut.sim1<- foo.complete.cut
foo.complete.cut.sim1$no.dem.branch<- 1
foo.complete.cut.sim1$second.session<- 0

foo.complete.cut.sim3<- foo.complete.cut
foo.complete.cut.sim3$no.dem.branch<- 3
foo.complete.cut.sim3$second.session<- 0

##include RE as known things because not predicting out of sample
##if forecast new year or value, put mean (gives more noise)
pred1<-  predict(mod.no.dem2b.lmer.cut, newdata = foo.complete.cut.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.cut, newdata = foo.complete.cut.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (-0.1238163)

###############
##If all moved to second session
foo.complete.cut.sim1<- foo.complete.cut
foo.complete.cut.sim1$no.dem.branch<- 1
foo.complete.cut.sim1$second.session<- 1

foo.complete.cut.sim3<- foo.complete.cut
foo.complete.cut.sim3$no.dem.branch<- 3
foo.complete.cut.sim3$second.session<- 1

##include RE as known things because not predicting out of sample
##if forecast new year or value, put mean (gives more noise)
pred1<-  predict(mod.no.dem2b.lmer.cut, newdata = foo.complete.cut.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.cut, newdata = foo.complete.cut.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (0.05422692)

###############################################
##Big.Cut2
foo.big.cut2<- data.frame(cbind(pct.change2$big.cut2, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.big.cut2)<- c("big.cut2", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.big.cut2<- na.omit(foo.big.cut2)

###############
##If all moved to first session
##most owned
foo.complete.big.cut2.sim1<- foo.complete.big.cut2
foo.complete.big.cut2.sim1$no.dem.branch<- 1
foo.complete.big.cut2.sim1$second.session<- 0

foo.complete.big.cut2.sim3<- foo.complete.big.cut2
foo.complete.big.cut2.sim3$no.dem.branch<- 3
foo.complete.big.cut2.sim3$second.session<- 0

##include RE as known things because not predicting out of sample
##if forecast new year or value, put mean (gives more noise)
pred1<-  predict(mod.no.dem2b.lmer.big.cut2, newdata = foo.complete.big.cut2.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.big.cut2, newdata = foo.complete.big.cut2.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (0.04867565)

###############
##If all moved to second session
##most owned
foo.complete.big.cut2.sim1<- foo.complete.big.cut2
foo.complete.big.cut2.sim1$no.dem.branch<- 1
foo.complete.big.cut2.sim1$second.session<- 1

foo.complete.big.cut2.sim3<- foo.complete.big.cut2
foo.complete.big.cut2.sim3$no.dem.branch<- 3
foo.complete.big.cut2.sim3$second.session<- 1

##include RE as known things because not predicting out of sample
##if forecast new year or value, put mean (gives more noise)
pred1<-  predict(mod.no.dem2b.lmer.big.cut2, newdata = foo.complete.big.cut2.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.big.cut2, newdata = foo.complete.big.cut2.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (-0.04593928)

#################################
remove(mod.no.dem2b.lmer.mess.with)
remove(mod.no.dem2b.lmer.cut)
remove(mod.no.dem2b.lmer.big.cut)
remove(mod.no.dem2b.lmer.big.cut2)
remove(mod.no.dem2b.lmer.big.increase)
remove(mod.no.dem2b.lmer.big.increase2)


#########################################
##  No Dem Branches                     #
##  Most owned                          #
##  First and second session interaction#
##  Table 2                             #
#########################################
##Dem most owned
##interaction with second session
##Leave alone or change it
mod.no.dem2b.lmer.owned.mess.with<- glmer(mess.with ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.dem2b.lmer.owned.cut<- glmer(cut ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.cut)

##If cut, small % cut or big out (1 is big cut)
##50% as big cut
mod.no.dem2b.lmer.owned.big.cut2<- glmer(big.cut2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.cut2)

##If increase, small % or big increase (1 if big)
##50% as big increase
mod.no.dem2b.lmer.owned.big.increase2<- glmer(big.increase2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.increase2)

###############################################
##Use forge version to rpedict move from 1 Dem branch to 3 Dem branches
##copy complete data for regression
##Cut
foo.cut<- data.frame(cbind(pct.change2$cut, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.cut)<- c("cut", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.cut<- na.omit(foo.cut)
##For Most owned
foo.complete.mo.cut<- subset(foo.complete.cut, foo.complete.cut$most.owned==1)

##simulations with all other variables at observed but change no.dem.branch from 1 to 3
##most owned
##If all moved to first session
foo.complete.mo.cut.sim1<- foo.complete.mo.cut
foo.complete.mo.cut.sim1$no.dem.branch<- 1
foo.complete.mo.cut.sim1$second.session<- 0

foo.complete.mo.cut.sim3<- foo.complete.mo.cut
foo.complete.mo.cut.sim3$no.dem.branch<- 3
foo.complete.mo.cut.sim3$second.session<- 0

pred1<-  predict(mod.no.dem2b.lmer.owned.cut, newdata = foo.complete.mo.cut.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.owned.cut, newdata = foo.complete.mo.cut.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (-0.2416)

###############
##If all moved to second session
##most owned
foo.complete.mo.cut.sim1<- foo.complete.mo.cut
foo.complete.mo.cut.sim1$no.dem.branch<- 1
foo.complete.mo.cut.sim1$second.session<- 1

foo.complete.mo.cut.sim3<- foo.complete.mo.cut
foo.complete.mo.cut.sim3$no.dem.branch<- 3
foo.complete.mo.cut.sim3$second.session<- 1

##include RE as known things because not predicting out of sample
pred1<-  predict(mod.no.dem2b.lmer.owned.cut, newdata = foo.complete.mo.cut.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.owned.cut, newdata = foo.complete.mo.cut.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (0.05614516)

###############################################
##Big.Cut2
foo.big.cut2<- data.frame(cbind(pct.change2$big.cut2, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.big.cut2)<- c("big.cut2", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.big.cut2<- na.omit(foo.big.cut2)
##For Most owned
foo.complete.mo.big.cut2<- subset(foo.complete.big.cut2, foo.complete.big.cut2$most.owned==1)

###############
##If all moved to first session
##most owned
foo.complete.mo.big.cut2.sim1<- foo.complete.mo.big.cut2
foo.complete.mo.big.cut2.sim1$no.dem.branch<- 1
foo.complete.mo.big.cut2.sim1$second.session<- 0

foo.complete.mo.big.cut2.sim3<- foo.complete.mo.big.cut2
foo.complete.mo.big.cut2.sim3$no.dem.branch<- 3
foo.complete.mo.big.cut2.sim3$second.session<- 0

##include RE as known things because not predicting out of sample
pred1<-  predict(mod.no.dem2b.lmer.owned.big.cut2, newdata = foo.complete.mo.big.cut2.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.owned.big.cut2, newdata = foo.complete.mo.big.cut2.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (0.1314378)

###############
##If all moved to second session
##most owned
foo.complete.mo.big.cut2.sim1<- foo.complete.mo.big.cut2
foo.complete.mo.big.cut2.sim1$no.dem.branch<- 1
foo.complete.mo.big.cut2.sim1$second.session<- 1

foo.complete.mo.big.cut2.sim3<- foo.complete.mo.big.cut2
foo.complete.mo.big.cut2.sim3$no.dem.branch<- 3
foo.complete.mo.big.cut2.sim3$second.session<- 1

##include RE as known things because not predicting out of sample
pred1<-  predict(mod.no.dem2b.lmer.owned.big.cut2, newdata = foo.complete.mo.big.cut2.sim1,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)

pred3<-  predict(mod.no.dem2b.lmer.owned.big.cut2, newdata = foo.complete.mo.big.cut2.sim3,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = TRUE)

##observed avergage treatment effect
mean(pred3 - pred1) ##on scale of DV (-0.07056066)


###################################
remove(mod.no.dem2b.lmer.owned.mess.with)
remove(mod.no.dem2b.lmer.owned.cut)
remove(mod.no.dem2b.lmer.owned.big.cut2)
remove(mod.no.dem2b.lmer.owned.big.increase2)


################################
## No Dem Branches Yearly     ##
## Dem most owned             ##
##Table 3                     ##
################################
##Regressions on yearly data
##Leave alone or change it
y.no.dem.2.mess.with<- lm(mess.with.pct.mo ~ no.dem.branch + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag  +
                      no.dem.branch*second.session,
                    data=yearly)
summary(y.no.dem.2.mess.with)

##If change it, cut or increase it (1 is cut)
y.no.dem.2.cut<- lm(cut.pct.mo ~ no.dem.branch + net.dem +
                 paygo + grh + unemploy + total_surplus_pct_lag +
                 no.dem.branch*second.session,
                    data=yearly)
summary(y.no.dem.2.cut)

##If cut, small % cut or zero out (1 is zero out)
y.no.dem.2.big.cut2<- lm(big.cut2.pct.mo ~ no.dem.branch + net.dem +
                     paygo + grh + unemploy + total_surplus_pct_lag +
                     no.dem.branch*second.session,
                    data=yearly)
summary(y.no.dem.2.big.cut2)

##If increase, small % or new (1 if new)
y.no.dem.2.big.increase2<- lm(big.increase2.pct.mo ~ no.dem.branch + net.dem +
                          paygo + grh + unemploy + total_surplus_pct_lag +
                    no.dem.branch*second.session,
                    data=yearly)
summary(y.no.dem.2.big.increase2)


#####################################################
##Robustness checks of models for appendix         ##
#####################################################
##Different cutpoints for BIG CUTS and BIG INCREASES#
#####################################################
##Define cutpoints 24, 40, 60, 75
##alternate definition of big change (baseline is greater than 50%)
pct.change2$big.increase25<- ifelse(pct.change2$new.event==1 |
                                   pct.change2$pctchangenom > 25.0, 1,
                             ifelse(pct.change2$pos.change.event==1, 0, NA))
pct.change2$big.cut25<- ifelse(pct.change2$zeroed.event==1 |
                                   pct.change2$pctchangenom < -25.0, 1,
                             ifelse(pct.change2$neg.change.event==1, 0, NA))

pct.change2$big.increase40<- ifelse(pct.change2$new.event==1 |
                                   pct.change2$pctchangenom > 40.0, 1,
                             ifelse(pct.change2$pos.change.event==1, 0, NA))
pct.change2$big.cut40<- ifelse(pct.change2$zeroed.event==1 |
                                   pct.change2$pctchangenom < -40.0, 1,
                             ifelse(pct.change2$neg.change.event==1, 0, NA))

pct.change2$big.increase60<- ifelse(pct.change2$new.event==1 |
                                   pct.change2$pctchangenom > 60.0, 1,
                             ifelse(pct.change2$pos.change.event==1, 0, NA))
pct.change2$big.cut60<- ifelse(pct.change2$zeroed.event==1 |
                                   pct.change2$pctchangenom < -60.0, 1,
                             ifelse(pct.change2$neg.change.event==1, 0, NA))

pct.change2$big.increase75<- ifelse(pct.change2$new.event==1 |
                                   pct.change2$pctchangenom > 75.0, 1,
                             ifelse(pct.change2$pos.change.event==1, 0, NA))
pct.change2$big.cut75<- ifelse(pct.change2$zeroed.event==1 |
                                   pct.change2$pctchangenom < -75.0, 1,
                             ifelse(pct.change2$neg.change.event==1, 0, NA))
##big.cut2 is 50 percent

###############################
##LMER models                ##
###############################
####################################
##No dem branches                 ##
##With all subaccounts            ##
##interaction with second session ##
##Table A1                        ##
####################################
##If cut, small % cut or big cut/zero out (1 is big cut)
##alternate definition
mod.no.dem2b.lmer.big.cut25<- glmer(big.cut25 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.cut25)

##If increase, small % or big/new (1 if big increase)
##alternate definition of big change
mod.no.dem2b.lmer.big.increase25<- glmer(big.increase25 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.increase25)

###################
##40%
##If cut, small % cut or big/zero out
mod.no.dem2b.lmer.big.cut40<- glmer(big.cut40 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.cut40)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem2b.lmer.big.increase40<- glmer(big.increase40 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.increase40)

###################
##60%
##If cut, small % cut or big/zero out
##alternate definition
mod.no.dem2b.lmer.big.cut60<- glmer(big.cut60 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.cut60)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem2b.lmer.big.increase60<- glmer(big.increase60 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.increase60)

########################
##75%
##If cut, small % cut or big/zero out
##alternate definition
mod.no.dem2b.lmer.big.cut75<- glmer(big.cut75 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.cut75)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem2b.lmer.big.increase75<- glmer(big.increase75 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.big.increase75)


remove(mod.no.dem2b.lmer.big.cut25)
remove(mod.no.dem2b.lmer.big.increase25)
remove(mod.no.dem2b.lmer.big.cut40)
remove(mod.no.dem2b.lmer.big.increase40)
remove(mod.no.dem2b.lmer.big.cut60)
remove(mod.no.dem2b.lmer.big.increase60)
remove(mod.no.dem2b.lmer.big.cut75)
remove(mod.no.dem2b.lmer.big.increase75)

#########################################
##  No Dem Branches                     #
##  Most owned                          #
##  First and second session interaction#
##  Table A2                            #
#########################################
##25%
##big cut
##alternate definition
mod.no.dem2b.lmer.owned.big.cut25<- glmer(big.cut25 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.cut25)

##If increase, small% or big/new
##alternate definition of big change
mod.no.dem2b.lmer.owned.big.increase25<- glmer(big.increase25 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.increase25)

########################
##40%
##big cut
##alternate definition
mod.no.dem2b.lmer.owned.big.cut40<- glmer(big.cut40 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.cut40)

##If increase, small% of big/new
##alternate definition of big change
mod.no.dem2b.lmer.owned.big.increase40<- glmer(big.increase40 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.increase40)

#########################
##60%
##big cut
##alternate definition
mod.no.dem2b.lmer.owned.big.cut60<- glmer(big.cut60 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.cut60)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem2b.lmer.owned.big.increase60<- glmer(big.increase60 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.increase60)

#################################
##75%
##big cut
##alternate definition
mod.no.dem2b.lmer.owned.big.cut75<- glmer(big.cut75 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.cut75)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem2b.lmer.owned.big.increase75<- glmer(big.increase75 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.lmer.owned.big.increase75)

remove(mod.no.dem2b.lmer.owned.big.cut25)
remove(mod.no.dem2b.lmer.owned.big.increase25)
remove(mod.no.dem2b.lmer.owned.big.cut40)
remove(mod.no.dem2b.lmer.owned.big.increase40)
remove(mod.no.dem2b.lmer.owned.big.cut60)
remove(mod.no.dem2b.lmer.owned.big.increase60)
remove(mod.no.dem2b.lmer.owned.big.cut75)
remove(mod.no.dem2b.lmer.owned.big.increase75)



###############################
##Yearly models              ##
##Tables A3 and A4           ##
###############################
##Most owned aggregated by year
##25% as big
yearly$big.cut25.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.cut25.pct.mo[i]<- (sum(pct.change2$big.cut25[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.cut25[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.cut25.pct.mo)

yearly$big.increase25.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.increase25.pct.mo[i]<- (sum(pct.change2$big.increase25[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.increase25[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.increase25.pct.mo)

##Regressions on yearly data
##alternate definition
y.no.dem.2.big.cut25<- lm(big.cut25.pct.mo ~ no.dem.branch + net.dem +
                     paygo + grh + unemploy + total_surplus_pct_lag +
                     no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.cut25)

##If increase, small % or big/new
##alternate definition of big change
y.no.dem.2.big.increase25<- lm(big.increase25.pct.mo ~ no.dem.branch + net.dem +
                          paygo + grh + unemploy + total_surplus_pct_lag +
                    no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.increase25)

############################################
##40% as big
yearly$big.cut40.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.cut40.pct.mo[i]<- (sum(pct.change2$big.cut40[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.cut40[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.cut40.pct.mo)

yearly$big.increase40.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.increase40.pct.mo[i]<- (sum(pct.change2$big.increase40[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.increase40[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.increase40.pct.mo)

##Regressions on yearly data
##alternate definition
y.no.dem.2.big.cut40<- lm(big.cut40.pct.mo ~ no.dem.branch + net.dem +
                     paygo + grh + unemploy + total_surplus_pct_lag +
                     no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.cut40)

##If increase, small % or big/new
##alternate definition of big change
y.no.dem.2.big.increase40<- lm(big.increase40.pct.mo ~ no.dem.branch + net.dem +
                          paygo + grh + unemploy + total_surplus_pct_lag +
                    no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.increase40)

##############################################
##60% as big
yearly$big.cut60.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.cut60.pct.mo[i]<- (sum(pct.change2$big.cut60[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.cut60[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.cut60.pct.mo)

yearly$big.increase60.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.increase60.pct.mo[i]<- (sum(pct.change2$big.increase60[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.increase60[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.increase60.pct.mo)

##Regressions on yearly data
##alternate definition
y.no.dem.2.big.cut60<- lm(big.cut60.pct.mo ~ no.dem.branch + net.dem +
                     paygo + grh + unemploy + total_surplus_pct_lag +
                     no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.cut60)

##If increase, small % or big/new
##alternate definition of big change
y.no.dem.2.big.increase60<- lm(big.increase60.pct.mo ~ no.dem.branch + net.dem +
                          paygo + grh + unemploy + total_surplus_pct_lag +
                    no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.increase60)

##########################################
##75% as big
yearly$big.cut75.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.cut75.pct.mo[i]<- (sum(pct.change2$big.cut75[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.cut75[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.cut75.pct.mo)

yearly$big.increase75.pct.mo<- NULL
for (i in 1:length(i.index)){
  yearly$big.increase75.pct.mo[i]<- (sum(pct.change2$big.increase75[pct.change2$year==i.index[i] & pct.change2$most.owned==1], na.rm=TRUE)/
                             length(na.omit(pct.change2$big.increase75[pct.change2$year==i.index[i] & pct.change2$most.owned==1])))*100
}
summary(yearly$big.increase75.pct.mo)

##Regressions on yearly data
##alternate definition
y.no.dem.2.big.cut75<- lm(big.cut75.pct.mo ~ no.dem.branch + net.dem +
                     paygo + grh + unemploy + total_surplus_pct_lag +
                     no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.cut75)

##If increase, small % or big/new
##alternate definition of big change
y.no.dem.2.big.increase75<- lm(big.increase75.pct.mo ~ no.dem.branch + net.dem +
                          paygo + grh + unemploy + total_surplus_pct_lag +
                    no.dem.branch*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2.big.increase75)


remove(y.no.dem.2.big.cut25)
remove(y.no.dem.2.big.increase25)
remove(y.no.dem.2.big.cut40)
remove(y.no.dem.2.big.increase40)
remove(y.no.dem.2.big.cut60)
remove(y.no.dem.2.big.increase60)
remove(y.no.dem.2.big.cut75)
remove(y.no.dem.2.big.increase75)


############################################
## Number Dem Branches as Factor          ##
############################################
##LMER models

####################################
##Number branches as factor       ##
##With all subaccounts            ##
##interaction with second session ##
##Table B1                        ##
####################################
##Leave alone or change it
mod.no.dem3b.lmer.mess.with<- glmer(mess.with ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.dem3b.lmer.cut<- glmer(cut ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.cut)

##If cut, small % cut or big/zero out
##alternate definition
mod.no.dem3b.lmer.big.cut2<- glmer(big.cut2 ~ as.factor(no.dem.branch)*second.session + net.dem +
                                  paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.big.cut2)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem3b.lmer.big.increase2<- glmer(big.increase2 ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.big.increase2)

remove(mod.no.dem3b.lmer.mess.with)
remove(mod.no.dem3b.lmer.cut)
remove(mod.no.dem3b.lmer.big.cut2)
remove(mod.no.dem3b.lmer.big.increase2)

#########################################
##  No Dem Branches as factor           #
##  Most owned                          #
##  First and second session interaction#
##  Table B2                            #
#########################################
##interaction with second session
##Leave alone or change it
mod.no.dem3b.lmer.owned.mess.with<- glmer(mess.with ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.owned.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.dem3b.lmer.owned.cut<- glmer(cut ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.owned.cut)

##If cut, small % cut or big/zero out
##alternate definition
mod.no.dem3b.lmer.owned.big.cut2<- glmer(big.cut2 ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.owned.big.cut2)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem3b.lmer.owned.big.increase2<- glmer(big.increase2 ~ as.factor(no.dem.branch)*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem3b.lmer.owned.big.increase2)

remove(mod.no.dem3b.lmer.owned.mess.with)
remove(mod.no.dem3b.lmer.owned.cut)
remove(mod.no.dem3b.lmer.owned.big.cut2)
remove(mod.no.dem3b.lmer.owned.big.increase2)

#################################
##Yearly model, most owned     ##
##Number Dem branches as factor##
##Table B3                     ##
#################################
##Regressions on yearly data
##Leave alone or change it
y.no.dem.2f.mess.with<- lm(mess.with.pct.mo ~ as.factor(no.dem.branch) + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag  +
                      as.factor(no.dem.branch)*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2f.mess.with)

##If change it, cut or increase it (1 is cut)
y.no.dem.2f.cut<- lm(cut.pct.mo ~ as.factor(no.dem.branch) + net.dem +
                 paygo + grh + unemploy + total_surplus_pct_lag +
                 as.factor(no.dem.branch)*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2f.cut)

##If cut, small % cut or big/zero out
##alternate definition
y.no.dem.2f.big.cut2<- lm(big.cut2.pct.mo ~ as.factor(no.dem.branch) + net.dem +
                     paygo + grh + unemploy + total_surplus_pct_lag +
                     as.factor(no.dem.branch)*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2f.big.cut2)

##If increase, small % or big/new
##alternate definition of big change
y.no.dem.2f.big.increase2<- lm(big.increase2.pct.mo ~ as.factor(no.dem.branch) + net.dem +
                          paygo + grh + unemploy + total_surplus_pct_lag +
                    as.factor(no.dem.branch)*second.session,
                    data=yearly,
                    #subset=(yearly$second.session==0)
                        )
summary(y.no.dem.2f.big.increase2)

remove(y.no.dem.2f.mess.with)
remove(y.no.dem.2f.cut)
remove(y.no.dem.2f.big.cut2)
remove(y.no.dem.2f.big.increase2)



#########################################
##  No Dem Branches continuous          #
##  Most owned as interaction           #
##  First and second session interaction#
##  Table C1                            #
#########################################
##interaction with second session
##Leave alone or change it
mod.no.interact.mo.dem4b.lmer.owned.mess.with<- glmer(mess.with ~ no.dem.branch*second.session*most.owned + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    family=binomial(link="logit"))
summary(mod.no.interact.mo.dem4b.lmer.owned.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.interact.mo.dem4b.lmer.owned.cut<- glmer(cut ~ no.dem.branch*second.session*most.owned + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    family=binomial(link="logit"))
summary(mod.no.interact.mo.dem4b.lmer.owned.cut)

##If cut, small % cut or big/zero out
##alternate definition
mod.no.interact.mo.dem4b.lmer.owned.big.cut2<- glmer(big.cut2 ~ no.dem.branch*second.session*most.owned + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    family=binomial(link="logit"))
summary(mod.no.interact.mo.dem4b.lmer.owned.big.cut2)


##If increase, small % or big/new
##alternate definition of big change
mod.no.interact.mo.dem4b.lmer.owned.big.increase2<- glmer(big.increase2 ~ no.dem.branch*second.session*most.owned + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    family=binomial(link="logit"))
summary(mod.no.interact.mo.dem4b.lmer.owned.big.increase2)


remove(mod.no.interact.mo.dem4b.lmer.owned.mess.with)
remove(mod.no.interact.mo.dem4b.lmer.owned.cut)
remove(mod.no.interact.mo.dem4b.lmer.owned.big.cut2)
remove(mod.no.interact.mo.dem4b.lmer.owned.big.increase2)


######################################
##  Robustness Checks               ##
##  Assess party control by branch  ##
######################################
##Note that writing over model names

###############################
##  All Subaccounts          ##
##  Drop no.dem.branch       ##
##  Use # Dem chambers + Dem ##
##  Pres                     ##
##  Table D1                 ##
###############################
##Drop no.dem.branch, but add dem.pres and no.dem.chambers
##Leave alone or change it
mod.no.dem2b.glmer.mess.with<- glmer(mess.with ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.glmer.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.dem2b.glmer.cut<- glmer(cut ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.glmer.cut)

##If cut, small % cut or big/zero out
##alternate definition
mod.no.dem2b.glmer.big.cut2<- glmer(big.cut2 ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.glmer.big.cut2)

##If increase, small % or big/new
##alternate definition of big change
mod.no.dem2b.glmer.big.increase2<- glmer(big.increase2 ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.no.dem2b.glmer.big.increase2)



###############################
## Dem Most Owned Subaccounts##
##  Drop no.dem.branch       ##
##  Use # Dem chambers + Dem ##
##  Pres                     ##
##  Table D2                 ##
###############################
##Drop no.dem.branch, but add dem.pres and no.dem.chambers
##Leave alone or change it
mod.mo.no.dem2b.glmer.mess.with<- glmer(mess.with ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.mo.no.dem2b.glmer.mess.with)

##If change it, cut or increase it (1 is cut)
mod.mo.no.dem2b.glmer.cut<- glmer(cut ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.mo.no.dem2b.glmer.cut)

##If cut, small % cut or big/zero out
##alternate definition
mod.mo.no.dem2b.glmer.big.cut2<- glmer(big.cut2 ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.mo.no.dem2b.glmer.big.cut2)

##If increase, small % or big/new
##alternate definition of big change
mod.mo.no.dem2b.glmer.big.increase2<- glmer(big.increase2 ~ no.dem.chamber*second.session + dem.pres*second.session + net.dem + dem.pres +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
summary(mod.mo.no.dem2b.glmer.big.increase2)



#####################################
##  Graph pct of cuts that are big ##
##  By variations on 2 dem branch  ##
##  Dem Owned Subaccounts          ##
##  Figure E1                      ##
#####################################
############
##Percent of cuts that are big
##Distribution of changes for Number Dem Branches

##Dem 1 branch
##Dem Pres
dem1p.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==1 & pct.change2$dem.pres==1 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==1 & pct.change2$dem.pres==1 & pct.change2$most.owned==1], na.rm=T)))*100
dem1p.pctbig.mo ##38.11

##Dem house
dem1h.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==1 & pct.change2$rep.house==0 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==1 & pct.change2$rep.house==0 & pct.change2$most.owned==1], na.rm=T)))*100
dem1h.pctbig.mo ##36.22

##Dem Senate
dem1s.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==1 & pct.change2$rep.senate==0 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==1 & pct.change2$rep.senate==0 & pct.change2$most.owned==1], na.rm=T)))*100
dem1s.pctbig.mo ##34.69

##2 Dem Branches
##House & Pres Dem (Senate Rep) NONE
dem2hp.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==2 & pct.change2$rep.senate==1 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==2 & pct.change2$rep.senate==1 & pct.change2$most.owned==1], na.rm=T)))*100
dem2hp.pctbig.mo

##Senate & Pres Dem (House Rep) NONE
dem2sp.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==2 & pct.change2$rep.house==1 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==2 & pct.change2$rep.house==1 & pct.change2$most.owned==1], na.rm=T)))*100
dem2sp.pctbig.mo

##House/Senate Dem + Rep Pres
dem2hs.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==2 & pct.change2$dem.pres==0 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==2 & pct.change2$dem.pres==0 & pct.change2$most.owned==1], na.rm=T)))*100
dem2hs.pctbig.mo ##53.29

##3Dem branches
dem3.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==3 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==3 & pct.change2$most.owned==1], na.rm=T)))*100
dem3.pctbig.mo ##50.12

##Original
dem1.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==1 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==1 & pct.change2$most.owned==1], na.rm=T)))*100
dem1.pctbig.mo

dem2.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==2 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==2 & pct.change2$most.owned==1], na.rm=T)))*100
dem2.pctbig.mo

dem3.pctbig.mo<- ((sum(pct.change2$big.cut2[pct.change2$no.dem.branch==3 & pct.change2$most.owned==1], na.rm=T)/sum(pct.change2$cut[pct.change2$no.dem.branch==3 & pct.change2$most.owned==1], na.rm=T)))*100
dem3.pctbig.mo


##barplot
par(mfrow=c(1,1), las=1)
barplot(c(dem1.pctbig.mo, dem1p.pctbig.mo, dem1h.pctbig.mo, dem1s.pctbig.mo,
          dem2.pctbig.mo, #dem2hp.pctbig.mo, dem2sp.pctbig.mo, dem2hs.pctbig.mo,
          dem3.pctbig.mo),
         ylim=c(0,70),
     main="", #Distribution of Budgetary Cuts
     xlab="Number of Democratic Institutions",
        ylab="Percent of Big Cuts",
     names.arg=c("1 (any)", "1 (Pres)", "1 (House)", "1 (Senate)",
     "2", #"2 Dem (HP)", "2 Dem (SP)", "2 Dem (HS)",
     "3"),
        cex.names=.7)


######################################################################
##  Bootstrap Main Models for Average Treatment Effect of Moving    ##
##  from 1 to 3 Democratically controlled institutions              ##
######################################################################
#########################################
##  No Dem Branches                     #
##  All                                 #
##  First and second session interaction#
#########################################
##Leave alone or change it
mod.no.dem2b.lmer.mess.with<- glmer(mess.with ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.dem2b.lmer.cut<- glmer(cut ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.cut)

##If cut, small % cut or zero out (1 is zero out)
##alternate definition
mod.no.dem2b.lmer.big.cut2<- glmer(big.cut2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.big.cut2)

##If increase, small % or new (1 if new)
##alternate definition of big change
mod.no.dem2b.lmer.big.increase2<- glmer(big.increase2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    #subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.big.increase2)


###############################################
##Use forge version to predict move from 1 Dem branch to 3 Dem branches
##copy complete data for bootstrap
##mess.with
foo.mess.with<- data.frame(cbind(pct.change2$mess.with, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.mess.with)<- c("mess.with", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.mess.with<- na.omit(foo.mess.with)


##cut
foo.cut<- data.frame(cbind(pct.change2$cut, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.cut)<- c("cut", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.cut<- na.omit(foo.cut)


##big.cut2
foo.big.cut2<- data.frame(cbind(pct.change2$big.cut2, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.big.cut2)<- c("big.cut2", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.big.cut2<- na.omit(foo.big.cut2)


##big.increase2
foo.big.increase2<- data.frame(cbind(pct.change2$big.increase2, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.big.increase2)<- c("big.increase2", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.big.increase2<- na.omit(foo.big.increase2)



##simulations with all other variables at observed but change no.dem.branch from 1 to 3

#################################
##        Mess.with             ##
#################################
##Average treatment effect of moving from 1 to 3 branches for all observed values (first session, then second session)
##bootstap
temp.mess.with<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.mess.with), replace=TRUE)
    sim.sample<- foo.complete.mess.with[sim.sample,] ##take rows from sample
    sim.model<- glmer(mess.with ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.mess.with[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.mess.with[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.mess.with, "glmer_bootstrap_messwith_all_13.csv", sep=",", row.names=FALSE)


#################################
##        cut                  ##
#################################
##Average treatment effect of moving from 1 to 2 branches for all observed values (first session, then second session)
##bootstap
temp.cut<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.cut), replace=TRUE)
    sim.sample<- foo.complete.cut[sim.sample,] ##take rows from sample
    sim.model<- glmer(cut ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.cut[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.cut[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.cut, "glmer_bootstrap_cut_all_13.csv", sep=",", row.names=FALSE)


#################################
##        big.cut2             ##
#################################
##Average treatment effect of moving from 1 to 2 branches for all observed values (first session, then second session)
##bootstap
temp.big.cut2<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.big.cut2), replace=TRUE)
    sim.sample<- foo.complete.big.cut2[sim.sample,] ##take rows from sample
    sim.model<- glmer(big.cut2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.big.cut2[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.big.cut2[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.big.cut2, "glmer_bootstrap_bigcut2_all_13.csv", sep=",", row.names=FALSE)


#################################
##     big.increase2           ##
#################################
##Average treatment effect of moving from 1 to 2 branches for all observed values (first session, then second session)
##bootstap
temp.big.increase2<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.big.increase2), replace=TRUE)
    sim.sample<- foo.complete.big.increase2[sim.sample,] ##take rows from sample
    sim.model<- glmer(big.increase2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.big.increase2[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.big.increase2[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.big.increase2, "glmer_bootstrap_bigincrease2_all_13.csv", sep=",", row.names=FALSE)


#########################################
##  No Dem Branches                     #
##  Most owned                          #
##  First and second session interaction#
#########################################
##Leave alone or change it
mod.no.dem2b.lmer.owned.mess.with<- glmer(mess.with ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.owned.mess.with)

##If change it, cut or increase it (1 is cut)
mod.no.dem2b.lmer.owned.cut<- glmer(cut ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.owned.cut)

##If cut, small % cut or zero out (1 is zero out)
##alternate definition
mod.no.dem2b.lmer.owned.big.cut2<- glmer(big.cut2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.owned.big.cut2)

##If increase, small % or new (1 if new)
##alternate definition of big change
mod.no.dem2b.lmer.owned.big.increase2<- glmer(big.increase2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=pct.change2,
                    subset=(pct.change2$most.owned==1),
                    family=binomial(link="logit"))
#summary(mod.no.dem2b.lmer.owned.big.increase2)


###############################################
##Use forge version to predict move from 1 Dem branch to 3 Dem branches
##copy complete data for bootstrap
##mess.with
foo.mess.with<- data.frame(cbind(pct.change2$mess.with, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.mess.with)<- c("mess.with", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.mess.with<- na.omit(foo.mess.with)
##For Most owned
foo.complete.mo.mess.with<- subset(foo.complete.mess.with, foo.complete.mess.with$most.owned==1)

##cut
foo.cut<- data.frame(cbind(pct.change2$cut, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.cut)<- c("cut", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.cut<- na.omit(foo.cut)
##For Most owned
foo.complete.mo.cut<- subset(foo.complete.cut, foo.complete.cut$most.owned==1)

##big.cut2
foo.big.cut2<- data.frame(cbind(pct.change2$big.cut2, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.big.cut2)<- c("big.cut2", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.big.cut2<- na.omit(foo.big.cut2)
##For Most owned
foo.complete.mo.big.cut2<- subset(foo.complete.big.cut2, foo.complete.big.cut2$most.owned==1)

##big.increase2
foo.big.increase2<- data.frame(cbind(pct.change2$big.increase2, pct.change2$no.dem.branch, pct.change2$second.session, pct.change2$net.dem, pct.change2$paygo,
                           pct.change2$grh, pct.change2$unemploy, pct.change2$total_surplus_pct_lag, pct.change2$dayspastfy, pct.change2$nsa,
                           pct.change2$most.owned))
names(foo.big.increase2)<- c("big.increase2", "no.dem.branch", "second.session", "net.dem", "paygo",
                           "grh", "unemploy", "total_surplus_pct_lag", "dayspastfy", "nsa",
                           "most.owned")
foo.complete.big.increase2<- na.omit(foo.big.increase2)
##For Most owned
foo.complete.mo.big.increase2<- subset(foo.complete.big.increase2, foo.complete.big.increase2$most.owned==1)


##simulations with all other variables at observed but change no.dem.branch from 1 to 3
##most owned

#################################
##        Mess.with            ##
#################################
##Average treatment effect of moving from 1 to 3 branches for all observed values (first session, then second session)
##bootstap
temp.mess.with<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.mo.mess.with), replace=TRUE)
    sim.sample<- foo.complete.mo.mess.with[sim.sample,] ##take rows from sample
    sim.model<- glmer(mess.with ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    subset=(sim.sample$most.owned==1),
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.mess.with[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.mess.with[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.mess.with, "glmer_bootstrap_messwith_most_owned_13.csv", sep=",", row.names=FALSE)


#################################
##        cut                  ##
#################################
##Average treatment effect of moving from 1 to 2 branches for all observed values (first session, then second session)
##bootstap
temp.cut<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.mo.cut), replace=TRUE)
    sim.sample<- foo.complete.mo.cut[sim.sample,] ##take rows from sample
    sim.model<- glmer(cut ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    subset=(sim.sample$most.owned==1),
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.cut[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.cut[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.cut, "glmer_bootstrap_cut_most_owned_13.csv", sep=",", row.names=FALSE)


#################################
##        big.cut2             ##
#################################
##Average treatment effect of moving from 1 to 2 branches for all observed values (first session, then second session)
##bootstap
temp.big.cut2<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.mo.big.cut2), replace=TRUE)
    sim.sample<- foo.complete.mo.big.cut2[sim.sample,] ##take rows from sample
    sim.model<- glmer(big.cut2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    subset=(sim.sample$most.owned==1),
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.big.cut2[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.big.cut2[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.big.cut2, "glmer_bootstrap_bigcut2_most_owned_13.csv", sep=",", row.names=FALSE)


#################################
##     big.increase2           ##
#################################
##Average treatment effect of moving from 1 to 2 branches for all observed values (first session, then second session)
##bootstap
temp.big.increase2<- matrix(nrow=1000, ncol=2) ##, mode="numeric") ##for simulated average treatment effects
for (i in 1:1000){
    sim.sample<- sample(1:nrow(foo.complete.mo.big.increase2), replace=TRUE)
    sim.sample<- foo.complete.mo.big.increase2[sim.sample,] ##take rows from sample
    sim.model<- glmer(big.increase2 ~ no.dem.branch*second.session + net.dem +
                      paygo + grh + unemploy + total_surplus_pct_lag + dayspastfy + (1|nsa),
                    data=sim.sample,
                    subset=(sim.sample$most.owned==1),
                    family=binomial(link="logit"))
    ##1 branch and first session
    sim.sample.sim1.first<- sim.sample
    sim.sample.sim1.first$no.dem.branch<- 1
    sim.sample.sim1.first$second.session<- 0
    ##2 branches and first session
    sim.sample.sim3.first<- sim.sample
    sim.sample.sim3.first$no.dem.branch<- 3
    sim.sample.sim3.first$second.session<- 0
    pred1.first<-  predict(sim.model, newdata = sim.sample.sim1.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.first<-  predict(sim.model, newdata = sim.sample.sim3.first,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##1 branch and second session
    sim.sample.sim1.second<- sim.sample
    sim.sample.sim1.second$no.dem.branch<- 1
    sim.sample.sim1.second$second.session<- 1
    ##2 branches and second session
    sim.sample.sim3.second<- sim.sample
    sim.sample.sim3.second$no.dem.branch<- 3
    sim.sample.sim3.second$second.session<- 1
    pred1.second<-  predict(sim.model, newdata = sim.sample.sim1.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    pred3.second<-  predict(sim.model, newdata = sim.sample.sim3.second,
    REform = NULL, terms = NULL,
    type = "response", allow.new.levels = FALSE)
    ##Average treatment effect for value of first or second session
    temp.big.increase2[i,1]<- mean(pred3.first - pred1.first) ##average treatment effect on scale of DV for first session
    temp.big.increase2[i,2]<- mean(pred3.second - pred1.second) ##average treatment effect on scale of DV for second session
}

##Output
write.table(temp.big.increase2, "glmer_bootstrap_bigincrease2_most_owned_13.csv", sep=",", row.names=FALSE)

#################################
##  Graph Average and 90% CI   ##
#################################
#################################
##Graph effects for all and most owned
##Plot bars with CI
library(plotrix)

##Bootstap results of moving from 1 to 3 branches (linear specification)
##All
mess.with.all<- read.csv("glmer_bootstrap_messwith_all_13.csv")
cut.all<- read.csv("glmer_bootstrap_cut_all_13.csv")
big.cut2.all<- read.csv("glmer_bootstrap_bigcut2_all_13.csv")
big.increase2.all<- read.csv("glmer_bootstrap_bigincrease2_all_13.csv")

##Most Owned
mess.with.mo<- read.csv("glmer_bootstrap_messwith_most_owned_13.csv")
cut.mo<- read.csv("glmer_bootstrap_cut_most_owned_13.csv")
big.cut2.mo<- read.csv("glmer_bootstrap_bigcut2_most_owned_13.csv")
big.increase2.mo<- read.csv("glmer_bootstrap_bigincrease2_most_owned_13.csv")

##Note: Need to multiple cut by -1 to get increase
increase.all<- cut.all*(-1)
increase.mo<- cut.mo*(-1)

##############
##Both (column 1 is first session effect, column 2 is second session effect
##90% CI All
##Figure 2
par(mfrow=c(2,1), las=1)
plotCI(1:4,
       c(mean(mess.with.all[,1]), mean(increase.all[,1]), mean(big.cut2.all[,1]), mean(big.increase2.all[,1])),
       ylim=c(-0.15,0.3),
       xlim=c(.5,4.5),
       xaxt="n",
       pch=1,
       main="(a) All Subaccounts",
       col="black",
       xlab="",
       ylab="Effect of of Dem. Control",
       cex.lab=.8,
       ui=c(quantile(mess.with.all[,1], 0.95), quantile(increase.all[,1], 0.95), quantile(big.cut2.all[,1], 0.95), quantile(big.increase2.all[,1], 0.95)),
       li=c(quantile(mess.with.all[,1], 0.05), quantile(increase.all[,1], 0.05), quantile(big.cut2.all[,1], 0.05), quantile(big.increase2.all[,1], 0.05)),
       err="y")
axis(1, at=seq(1,4,by=1), labels=c("Change", "Increase", "Big Cut", "Big Increase"), cex.axis=.8)

plotCI(c(1.1,2.1,3.1,4.1),
       c(mean(mess.with.all[,2]), mean(increase.all[,2]), mean(big.cut2.all[,2]), mean(big.increase2.all[,2])),
       col="grey",
       add=T,
       ui=c(quantile(mess.with.all[,2], 0.95), quantile(increase.all[,2], 0.95), quantile(big.cut2.all[,2], 0.95), quantile(big.increase2.all[,2], 0.95)),
       li=c(quantile(mess.with.all[,2], 0.05), quantile(increase.all[,2], 0.05), quantile(big.cut2.all[,2], 0.05), quantile(big.increase2.all[,2], 0.05)),
       err="y")

abline(h=0, lty=2, add=T)

legend(.5, -0.05,c ("First Session", "Second Session"), lty=c(1,1), col=c("black", "grey"), cex=.8, bty="n")


##Most owned subaccounts
##90% CI MO
plotCI(1:4,
       c(mean(mess.with.mo[,1]), mean(increase.mo[,1]), mean(big.cut2.mo[,1]), mean(big.increase2.mo[,1])),
       ylim=c(-0.15,0.3),
       xlim=c(.5,4.5),
       xaxt="n",
       pch=1,
       main="(b) Democratically Owned Subaccounts",
       col="black",
       xlab="",
       ylab="Effect of Dem. Control",
       cex.lab=.8,
       ui=c(quantile(mess.with.mo[,1], 0.95), quantile(increase.mo[,1], 0.95), quantile(big.cut2.mo[,1], 0.95), quantile(big.increase2.mo[,1], 0.95)),
       li=c(quantile(mess.with.mo[,1], 0.05), quantile(increase.mo[,1], 0.05), quantile(big.cut2.mo[,1], 0.05), quantile(big.increase2.mo[,1], 0.05)),
       err="y")
axis(1, at=seq(1,4,by=1), labels=c("Change", "Increase", "Big Cut", "Big Increase"), cex.axis=.8)

plotCI(c(1.1, 2.1,3.1,4.1),
       c(mean(mess.with.mo[,2]), mean(increase.mo[,2]), mean(big.cut2.mo[,2]), mean(big.increase2.mo[,2])),
       col="grey",
       add=T,
       ui=c(quantile(mess.with.mo[,2], 0.95), quantile(increase.mo[,2], 0.95), quantile(big.cut2.mo[,2], 0.95), quantile(big.increase2.mo[,2], 0.95)),
       li=c(quantile(mess.with.mo[,2], 0.05), quantile(increase.mo[,2], 0.05), quantile(big.cut2.mo[,2], 0.05), quantile(big.increase2.mo[,2], 0.05)),
       err="y")

abline(h=0, lty=2, add=T)

legend(.5, -0.05,c ("First Session", "Second Session"), lty=c(1,1), col=c("black", "grey"), cex=.8, bty="n")














